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ABSTRACT 

We study how the structure and variabiHty of magnetohydrodynamic (MHD) turbu- 
lence in accretion discs converge with domain size. Our results are based on a series of 
vertically stratified local simulations, computed using the Athena MHD code, that have 
fixed spatial resolution, but varying radial and azimuthal extent (from AR = 0.5H to 
16H, where H is the vertical scale height). We show that elementary local diagnostics 
of the turbulence, including the Shakura-Sunyaev a parameter, the ratio of Maxwell 
stress to magnetic energy, and the ratio of magnetic to fluid stresses, converge to 
within the precision of our measurements for spatial domains of radial size L^ > 2H. 
We obtain a ~ 0.02 — 0.03, consistent with other recent determinations. Very small 
domains {L^ = 0.5H) return anomalous results, independent of spatial resolution. 
This convergence with domain size, however, is only valid for a limited set of diagnos- 
tics: larger spatial domains admit the emergence of dynamically important mesoscale 
structures. In our largest simulations, the Maxwell stress shows a significant large 
scale non-local component, while the density develops long-lived axis5Tnmetric per- 
turbations ("zonal flows") at the 20% level. Most strikingly, the variability of the disc in 
fixed-sized patches decreases strongly as the simulation volume increases, while vari- 
ability in the magnetically dominated corona remains constant. Comparing our largest 
local simulations to global simulations with comparable spatial resolution, we find 
generally good agreement. There is no direct evidence that the presence of curvature 
terms or radial gradients in global calculations materially affect the turbulence, except 
to perhaps introduce an outer radial scale for mesoscale structures. The demonstrated 
importance of mean magnetic fields - seen in both large local and global simulations 
- implies, however, that the growth and saturation of these fields is likely of critical 
importance for the evolution of accretion discs. 
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1 INTRODUCTION 

Accretion of gas through a disc onto a central object plays a 
pivotal role in the formation of protostars and planets, and 
furnishes the only probe of the environment around black 
holes to date. If the disc is sufficiently ionized, accretion is 
driven by a powerful instability that is present in orbiting, 
magnetized gas: the magnetorotational instability (MRI; Bal- 
bus & Hawley 1998, 1991; Blaes & Balbus 1994; Goodman 
& Xu 1994). The MRI gives rise to sustained magnetohydro- 
djmamic (MHD) turbulence in which the turbulent fluctua- 
tions are correlated to produce a positive stress tensor and 
thus the outward angular momentum transport necessary for 
the gas to accrete. While numerous analytic investigations of 



the MRI's linear regime (e.g., Balbus & Hawley 1991; Blaes & 
Balbus 1994) have been insightful, a complete understanding 
of turbulent disc accretion requires the ability to probe the 
fully nonlinear behavior of this instability. This is best done 
through the use of numerical simulations (though, some ana- 
lytic work in the non-linear regime has been done; Goodman 
& Xu 1994; Pessah & Goodman 2009). 

Numerical studies of accretion discs can be either global 
or local. These classes of simulations are distinguished by the 
different length scales that are captured. A local patch of ac- 
cretion disc can be characterized by two length scales: the 
distance from the central object, _Ro and the local scale height 
(disc thickness), H ^ Cs/il (where Cs and fl are the sound 
speed and angular velocity at Ro, respectively). If we asso- 
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ciate the domain of a simulation with a length scale, L, then 
local simulations are characterized by i? ~ L <C -Ro, whereas 
global simulations are characterized hy H <§:. L '^ Rq. 

Local simulations invariably employ the shearing box ap- 
proximation (e.g., Hawley et al. 1995; Brandenburg et al. 
1995), in which the MHD equations are solved in the frame 
of a local, co-rotating patch of the disc. The assumption that 
H ^ L <^ Ro allows the use of a Cartesian geometry while 
retaining the essential dynamics of differential rotation. This 
simplifies the problem down to its basic ingredients, orbital 
shear and magnetization, allowing one to answer basic ques- 
tions about MRI-driven turbulence, such as "what determines 
the amplitude of turbulent fluctuations?" At fixed computa- 
tional cost, local simulations provide the best resolution of 
small-scale features in the turbulence, which is particularly 
advantageous if one wishes to study processes such as resistiv- 
ity or viscosity that are physically important only for Z <; H. 
Local simulations can in principle be simplified yet further 
by ignoring the vertical component of gravity ("unstratified" 
simulations) . However, in our current work, we only consider 
stratified calculations which include the vertical gravity from 
the central object. 

Global simulations (e.g., Armitage 1998; Hawley 2000), 
on the other hand, evolve a domain that is comparable to 
the distance from the central object. These simulations are 
not only t5^ically larger in spatial extent than local models - 
which may be important if for example long wavelength ra- 
dial and azimuthal fluctuations are present in the disc - they 
also include new physical effects. In particular, global mod- 
els include terms associated with disc curvature, and possess 
the asymmetry between inward and outward radial directions 
that is needed to explicitly represent mass accretion. How- 
ever, these advantages come at the expense of non-trivial "in- 
ner" and "outer" boundary conditions. For example, with the 
exception of black hole disk simulations, where an outflow 
boundary is applied at the inner boundary to account for ma- 
terial plummeting into an event horizon, particular methods 
(e.g., viscous or magnetic damping) have to be employed to 
reduce spurious effects from these boundaries. Global simu- 
lations also suffer from poor resolution of small scales, such 
that most global simulations to-date fail to meet the minimum 
requirements for numerical convergence that have been de- 
rived from local models (Hawley et al. 2011). Only recently 
have increases in computational power and algorithmic im- 
provements made it possible to compute global simulations 
at the resolutions required to adequately resolve the MRI at 
the level done in local simulations (Noble et al. 2010; Beck- 
with et al. 2011; Sorathia et al. 2011). 

The near congruence of resolutions between local and 
global simulations makes possible a comparison of the turbu- 
lent properties in these two regimes. One question is whether 
a local patch of a well-resolved global simulation resembles 
an equivalently resolved local simulation? If not, are any dis- 
crepancies due solely to the different sizes of the domains, 
or are they instead attributable to truly global effects such 
as curvature? One approach to addressing this question is to 
inspect the properties of local sub-domains within a global 
simulation and compare the statistical properties of these sub- 
domains to results from shearing boxes (Sorathia et al. 2010, 
2011). These studies suggest that some of the basic properties 
of local MRI turbulence, such as the relationship between tur- 
bulent stress and magnetic flux through a local patch (Hawley 



et al. 1995) and the tilt angle of the magnetic field autocorre- 
lation function (Guan et al. 2009), carry over to global calcu- 
lations. 

A complementary approach is to examine the properties 
of turbulence in large "mesoscale" simulations, characterized 
hy H <^ L <^ Ro (Guan & Gammie 2011). These simulations, 
which continue to utilize the shearing box approximation and 
are in this sense still local, allow study of the influence of long 
wavelength azimuthal and radial modes on the properties 
of the turbulence, while also capturing small scale turbulent 
fluctuations. Compared to global simulations, the mesoscale 
neglects the effect of curvature terms on the disc, which, in 
principle, allows the influence of these terms to be elucidated. 
Guan & Gammie (2011) (see also Davis et al. 2010) find that, 
in the mesoscale regime, turbulence is characterized by small 
scale fluctuations regardless of domain size. These authors 
also find that, away from the mid-plane, in the coronal re- 
gion, the magnetic field is correlated on scales of ~ IQH but 
does not contribute significantly to angular momentum trans- 
port. Larger scales are not, however, unimportant. Nelson & 
Gressel (2010) (see also Gressel et al. 2011; Yang et al. 2011) 
demonstrate that long wavelength (again » H) density cor- 
relations determine the amplitude of stochastic density fluc- 
tuations (spiral density waves), while the calculations of Jo- 
hansen et al. (2009) show long lived density/pressure zonal 
flows in geostrophic balance. 

Our goal in this work is to systematically study how the 
properties of disc turbulence change as one transitions from 
the local to the mesoscale regime. We do so though a series 
of shearing box simulations of increasing domain size, from 
i^ = Q.^H X Ly = 2H X L:, = 8H to L^ = 16H x Ly = 
32H xLz — 8H (here x, y and z refer to the radial, azimuthal, 
and vertical directions respectively) ^ . In particular, we wish to 
answer the following questions. How do various properties of 
MRI-driven turbulence depend on box size? As we push to- 
wards the mesoscale regime, what if any phenomena emerge 
at these larger scales? Is there a scale at which the turbulent 
properties resemble those in a global disc calculation? 

The rest of this paper is structured as follows. In §2, we 
give details of the numerical algorithm used to integrate the 
equations of ideal MHD in the local limit and the initial and 
boundary conditions used for the simulations. In §3, we de- 
tail the properties of the quasi-stationary turbulent state via 
volume-averaged quantities and the vertical structure of the 
turbulence. In §4, we examine the properties of mesoscale 
structures within the magnetic field and the implications for 
the locality of angular momentum transport. The derived vari- 
ability of accretion is examined in §5. In §6 we compare a 
subset of our results to those obtained from well-resolved 
global simulations. Finally, in §7, we summarize our results 
and discuss implications for the physics of magnetized accre- 
tion discs. 



^ Whether the neglect of radial gradients and curvature terms in our 
large local boxes is physically justified depends upon the geometric 
thickness of the disc being modeled: a domain of radial scale ^ H 
remains "local" in thin regions of an AGN disc where H/R < 10~^, 
whereas the same would not be true of a thicker protoplanetary disc. 
Our primary interest in the large boxes is as well-controlled model 
systems, though the actual value of H/R should be borne in mind 
when comparing our results to global simulations. 
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2 METHOD 

2.1 Numerical Algorithm 

For these calculations, we use Athena, a second-order accu- 
rate Godunov flux-conservative code for solving the equations 
of MHD. Athena uses the dimensionally unsplit corner trans- 
port upwind (CTU) method of Colella (1990) coupled with 
the third-order in space piecewise parabolic method (PPM) of 
Colella & Woodward (1984) and a constrained transport (CT; 
Evans & Hawley 1988) algorithm for preserving the V B = 
constraint. We use the HLLD Riemann solver to calculate the 
numerical fluxes (Miyoshi & Kusano 2005; Mignone 2007). A 
detailed description of the Athena algorithm and the results of 
various test problems are given in Gardiner & Stone (2005), 
Gardiner & Stone (2008), and Stone et al. (2008). 

The simulations use the shearing box approximation, a 
model for a local co-rotating disc patch whose size is small 
compared to the radial distance from the central object, Rq. 
We construct a local Cartesian frame, x = {R — Ro), y = Ro(l>, 
and z, co-rotating with an angular velocity Q, corresponding 
to the orbital frequency at Ro, the center of the box. In this 
frame, the equations of motion become (Hawley et al. 1995) : 



Table 1 . Shearing Box Simulations 
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2qpQ X — pQ, z — 2fi X pv 
dB 



V X (v X B) = 



where p is the mass density, pv is the momentum density, B 
is the magnetic field, P is the gas pressure, and q is the shear 
parameter, defined as g = — dlnfi/dlni?. We use q = 3/2, 
appropriate for a Keplerian disc. We assume an isothermal 
equation of state P = pc^, where Cg is the isothermal sound 
speed. From left to right, the source terms in the momentum 
equation correspond to radial tidal forces (gravity and cen- 
trifugal), vertical gravity, and the Coriolis force. Note that our 
system of units has the magnetic permeability ^ = 1. 

The numerical integration of the shearing box equations 
require additions to the Athena algorithm, the details of which 
can be found in Stone & Gardiner (2010) and in the Appendix 
of Simon et al. (2011). Briefly, we utilize Crank-Nicholson dif- 
ferencing to conserve epicyclic motion exactly and orbital ad- 
vection to subtract off the background shear flow (Stone & 
Gardiner 2010). The y boundary conditions are strictly peri- 
odic, whereas the x boundaries are shearing periodic (Haw- 
ley et al. 1995; Simon et al. 2011). The vertical boundaries 
are the outflow boundary conditions described in Simon et al. 
(2011). 



2.2 Simulation Parameters and Initial Conditions 

All of our simulations are vertically stratified, with an initial 
density corresponding to isothermal hydrostatic equilibrium. 



p{^, y, z) = poexp 



m 



(2) 



Label 




Domain Size 




Resolution 


Initial Field 




{L./- 


H X Ly/H X L^/H) 


(zones/H) 


Configuration 


FT0.5 




0.5x 2x 8 




32 


Flux tube 


FTO.Sm 




0.5 X 2 X 8 




72 


Flux tube 


FTO.Sh 




0.5 X 2 X 8 




144 


Flux tube 


FT2 




2x4x8 




32 


Flux tube 


FT4 




4x8x8 




32 


Flux tube 


FT8 




8 X 16 X 8 




32 


Flux tube 


Y16 




16 X 32 X 8 




36 


Toroidal 



Here, the scale height is defined as H = v^Cs/Q. 

where po = 1 is the mid-plane density, and H is the scale 
height in the disc. 



H = 



V2c 



(3) 



The isothermal sound speed, Cs = 7.07 x 10""', corresponding 
to an initial value for the mid-plane gas pressure of Po = 5 x 
10"'^. With Q. = 0.001, the value for the scale height isH = 1. 
A density floor of 10""' (10^'' for the smallest domain run) is 
applied to the physical domain as too small a density leads to 
a large Alfven speed and a very small timestep. Furthermore, 
numerical errors make it difficult to evolve regions of very 
small plasma /3. 

The initial magnetic field configuration for most of our 
calculations is the twisted azimuthal flux tube of Hirose et al. 
(2006), with minor modifications to the dimensions and /3 
values. In particular, the initial toroidal field. By, is given by 



By 



2Po 

/3y 



(BI + BI) 







if B^ + Bf / 
if B^ -f- B^ = 



(4) 



with Py — 100. The poloidal field components, B^ and B^, 
are calculated from the y component of the vector potential. 



Ay = 



1 4- cos 







(S)] 



if r < Ro 
if r > Bo 



(5) 



where r = ^/x'^ + z'^ and /3p — 1600 is the poloidal field /3 
value. We choose Bo to always be one fourth of the radial 
domain size; Bo — Lx/4:. Figure 1 shows a rendering of the 
initial magnetic field configuration from FT4. 

For the largest domain size, the aspect ratio is such that 
Bo = I/z/2. As was evident from an initial run at this size 
and with the flux tube field geometry, the flux tube cannot 
be properly contained within the domain without initial and 
dominant boundary condition effects. Thus, we chose our 
largest domain run, Y16, to be initialized with a constant 
13 = 100 toroidal field throughout the entire domain. For all of 
our calculations, random perturbations are added to the den- 
sity and velocity components to seed the MRI. Table 1 lists 
the simulations that we have analyzed for this work. The la- 
bel for each calculation describes the initial field geometry 
and the x domain size. So, for run FT0.5, FT denotes the flux 
tube geometry, and the 0.5 denotes Lx — 0.5B. As shown in 
the table, we have also run two higher resolution simulations 
equivalent to FT0.5; FT0.5m is carried out at 72 zones per H 
and FT0.5h is carried out at 144 zones per H. 
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Figure 1. Rendering of tiie initial magnetic field lines for the FT4 
run. This is the twisted azimuthal flux tube of Hirose et al. (2006), 
and it is the initial field configuration for all of our runs except for 
Y16. 



3 CHARACTERIZING THE TURBULENT STEADY STATE 

Although the focus of the current study is on the convergence 
of the properties of disc turbulence with varying domain size, 
this can only be meaningful if the simulations being compared 
have adequate spatial resolution. To satisfy ourselves that this 
is the case, we first compute the quality parameter Qj (Haw- 
ley et al. 2011), which characterizes the effective numerical 
resolution of the turbulence (i.e., how well-resolved the tur- 
bulence is). For a direction j — y, z, the quality parameter is 
defined as. 



2-r\vk\ . 



nAa 



\va\ = 



Ml 



(6) 



where Axj is the grid cell spacing along direction j. Here, the 
brackets denote a volume average for all x and y and for jz] < 
2H, and the overbar is a time average from orbit 50 onwards. 
Sano et al. (2004) demonstrate that Q ~ 6 is required in 
order for the MRI to be well-resolved in the saturated state 
(though this value is likely to be code-dependent). The values 
of Q measured for our simulations are given in the final two 
columns of Table 2. All of the simulations are well-resolved 
according to this criterion. 

Turning to physical rather than numerical quantities, we 
examine the total (Maxwell plus Rejniolds) stress normalized 
to the gas pressure. The time evolution of this measure, which 
is equivalent to the a-parameter introduced by Shakura & 
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Figure 2. Time history of volume-averaged Maxwell and Re5molds 
stresses normalized by the volume-averaged pressure for each shear- 
ing box domain size. The volume average is done for all x and y and 
for \z\ < 2H. The saturation level is roughly consistent between all 
of the calculations. 



Table 2. Saturation Characteristics 



Run 



MR BR^ 



BR„ 



BR^ 



Wy ^ z 



FT0.5 


0.022 


0.16 


4.6 


0.056 


0.92 


0.022 


67 


10 


FTO.Sm 


0.030 


0.13 


4.7 


0.057 


0.92 


0.024 


188 


30 


FTO.Sh 


0.032 


0.22 


5.5 


0.096 


0.86 


0.045 


297 


67 


FT2 


0.027 


0.38 


3.8 


0.11 


0.84 


0.051 


44 


11 


FT4 


0.023 


0.40 


3.5 


0.11 


0.84 


0.048 


39 


9 


FT8 


0.020 


0.41 


3.4 


0.11 


0.85 


0.045 


35 


8 


Y16 


0.024 


0.40 


3.5 


0.12 


0.83 


0.050 


44 


11 



Here, a is the ratio of stress to gas pressure (i.e., the traditional 
Shakura & Sunyaev (1973) a parameter), cimag is the ratio of 
Maxwell stress to magnetic energy, MR is the Maxwell to Reynolds 
stress ratio, BRi is the ratio of magnetic energy component i to the 
total magnetic energy, and Qi is the quality parameter along direction 
i as defined via equation (6) . 



Sunyaev (1973), is shown in Figure 2, where we have volume- 
averaged the data over all x and y and for \z\ < 2H. Following 
an initial transient period, which lasts for the first ^ 50 orbits 
of the evolution, the normalized stress levels are all roughly 
equal between different domain sizes. There is, however, a 
dramatic change in temporal variability as one goes to larger 
box sizes. This effect was present in the simulations of Davis 
et al. (2010), and we will discuss it in detail in §5. 

The data of Figure 2 allow us to define a time-period, af- 
ter the 50-orbit mark, where the flow is in a quasi-stationary 
turbulent state. We use the simulation data from this period to 
construct various time- and volume-averaged measures of the 
turbulence (where the volume-average is again performed for 
all X and y and for |z;| < 2H). We consider the following diag- 
nostics that characterize the physical state of the turbulence: 
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{pVxdVy — B^By) 



-BxBy) 



MR; 



(£2) 



(pVxSVy) 



BR, 






(7) 



These are, in turn, the traditional Shakura & Sunyaev (1973) 
a parameter; the ratio of Maxwell stress to magnetic energy, 
cfmag (Hawley et al. 2011)^; the ratio of Maxwell to Reynolds 
stress, MR; and the ratio of magnetic energy component i to 
the total magnetic energy, BR;. The angled brackets denote 
a volume average, whereas the overbars indicate a time av- 
erage. These quantities are shown in Figure 3 and in Table 
2. We find that, within the precision possible given the lim- 
ited duration of our simulations, all of these metrics appear 
to converge for domain sizes Lj, > 2H, i.e. for all but the 
smallest domain. We find that a varies between 0.02 — 0.03; 
cfmag ~ 0.4, consistent with expectations for converged MRI 
turbulence (Hawley et al. 2011; Sorathia et al. 2011); the ra- 
tio of Maxwell to Reynolds stresses is ~ 3 — 4; while finally 
BR^ ~ 0.11, BRy ~ 0.84, BR^ ~ 0.05. These diagnostics 
can be directly compared with the results of Beckwith et al. 
(2011), who characterized the properties of MRI-driven MHD 
turbulence in a well-resolved global simulation. In that work, 
the authors found that disc-turbulence measured in regions 
well away from the innermost stable circular orbit was char- 
acterized by a ~ 0.025 (directly in the middle of the range 
found here), a^ag ^ 0.3 and finally BR^ ~ 0.1, BRj, ~ 0.88, 
BRz ~ 0.02. We conclude that the basic measures of the tur- 
bulence reported here are broadly consistent with those in 
Beckwith et al. (2011). 

While boxes larger than 2H x AH x 8H show high lev- 
els of convergence in these diagnostics, the behaviour of the 
smallest box (FT0.5) is quite different. While a = 0.022 in 
this case, consistent with the other domain sizes presented 
here, the other diagnostics are significantly outside of the 
range of parameter space occupied by the larger domains. 
Most tellingly, we find that Omag is a factor 2.5 smaller than 
the expected value for converged MRI turbulence (Hawley 
et al. 2011; Sorathia et al. 2011). We have verified this be- 
havior for calculations with the same domain size but with 
72 and 144 zones per H (see Table 2). The run FTO.Sh ap- 
pears to have values closer to the converged values. However, 
analyzing the time history of these quantities show that they 
are not quite saturated between orbit 50 and the end of the 
calculation. In particular, Qmag is still decreasing during the 
period over which we average. Given limited computational 
resources, we did not run this case out further, but instead 
we expect that Qmag will saturate at a value lower than that 
presented here. The bottom line is the anomalous behavior of 
the small domain simulations is not a resolution effect. 

The turbulence in the quasi-stationary state can be fur- 
ther characterized by the vertical structure of the simulated 
disc. In Simon et al. (2011), we examined the statistical dis- 
tribution of turbulent velocity fluctuations as a function of 
height above the midplane in runs FT2, FT4, FT8 and Y16, 
and showed that there was a non-negligible supersonic com- 
ponent for \z\ > 3H (see Fig. 3 of Simon et al. 2011). In run 



^ This is closely related to the "tilt-angle" measured by several au- 
thors (Guan et al. 2009; Beckwith et al. 2011; Sorathia et al. 2011). 
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Figure 3. Several time- and volume-averaged quantities versus radial 
domain size. The left panel displays a (squares), amag (asterisks), 
and MR (x's), whereas the right panel displays BRy (squares), BR^; 
(triangles), and BR^ (diamonds). These diagnostics are defined by 
equation (7) and converge for L^/H > 2. 



FT2 and FT4, this component occupied ~ 1% and ~ 5% of 
the probability distribution space respectively, while for sim- 
ulations FT8 and Y16, the contribution of supersonic veloc- 
ity fluctuations converged at ~ 10% of the probability distri- 
bution space, highlighting the possibility of observable differ- 
ences between local and mesoscale treatments of magnetized 
accretion disc turbulence. 

Fig. 4 shows the vertical structure of the total stress nor- 
malized to the initial mid-plane gas pressure and the gas 
/? — 2P/B^ parameter in the quasi-stationary state. Here, 
the numerator and denominator for these quantities are time- 
(from orbit 50 onward) and horizontally-averaged separately 
before computing their ratio. With the exception of the small- 
est domain run, the vertical structure appears reasonably con- 
sistent between the various domain sizes; the stress is rela- 
tively flat until \z\ ~ 2H, after which it drops off rapidly, 
and \z\ ~ 2H is also where the /? ~ 1 line is crossed. The 
smallest domain again appears pathological in comparison 
to the larger domain sizes: the stress is sharply concentrated 
towards the midplane, where the disc is weakly magnetized 
(/? ~ 20 cf /? ~ 40 in the larger domains), while /3 < 1 for 

Fig. 5 shows the vertical profile of the energy density 
in toroidal, radial and vertical magnetic fields normalized to 
the total magnetic energy density, e.g. (Sf ) / ((_B^)) (i = 
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100.0 



10.0 




Figure 4. Time- and horizontally-averaged vertical profiles of the to- 
tal stress normalized to the initial mid-plane gas pressure (top panel) 
and gas f} = 2PjB'^ parameter (bottom panel). In both cases, the nu- 
merator and denominator are separately averaged. The green curve 
is from FT0.5, black is FT2, red is FT4, blue is FT8, and orange is 
Y16. With the exception of the smallest domain, the vertical structure 
of these quantities are consistent between all calculations: the stress 
is relatively flat until \z\ ~ 2H, after which it drops off rapidly, and 
\z\ ~ 2H is also where the /3 ~ 1 line is crossed. 



x, y, z). As before, we average these quantities separately be- 
fore taking their ratio and we omit data from the smallest 
domain (run FT0.5) for the purposes of clarity. As was the 
case with the data shown in Figure 4, the vertical structure of 
these quantities appear to be reasonably consistent between 
the various domain sizes: the energy density in toroidal fields 
dominates throughout the domain. In the region \z\ < IH, 
the ratio of energy densities is 4 : 40 : \ {x : y : z), while in 
the region \z\ > "iH, the energy density in vertical fields ex- 
ceeds that in radial fields, becoming comparable to (but still 
smaller than) the energy density in toroidal fields. This leads 
us to an interesting point: the largest domain (Y16) has ver- 
tical magnetic field strengths at 1^1 = \H that exceed the 
strength of the radial magnetic field here by a factor ~ 2 
and are within a factor of ~ 2 of the toroidal magnetic field 
strength. Given the vertical gradients of each component of 
the magnetic field strength, it appears that the energy density 
in vertical fields will exceed that in toroidal magnetic fields 
by |2| = 6-ff; such a configuration could lead to angular mo- 
mentum transport loss through a magnetically launched wind 



A 
A 

tN 

QD 
V 
V 



QD 
V 



10.000 


; 


; 






- 










>^'-^r^'''' ' '"^^S^'. --v.. 




1.000 


/-'■' -^-rj^ 


- 




/' ^'^s. 


■ 










j if "V-\ 


■ 




/ .' "W 








■ 


0.100 


' 1 i>:: =„.,=::-^. 15 


- 




/*- //■■'■■■ ■S.V^, V 






■ // /^■<*. ■ ■>>■'. 






/■? //>:-■■ ■-. ■vv^\ 






/ * fi ■*■ ■:\''<i.. 






■i /// ■% 


*'\\^ ■ 




1 o"/ ->?; 


%'< 




//■■/?'' N'» 


':■■ \ 


0.010 


- '■' ,' \ 


l---h 






\\ '■- 




;/'/ 


- '?■-= 




^11 ' 


\ ^\. 




''/ 


\ M-., 




■// 


\ n 




[/ 


jT 


0.001 


' . , . 1 , . , 1 . , . 1 , . 






z (H) 



Figure 5. The time- and horizontally-averaged vertical profile of the 
magnetic field strength in individual components, (B^) (i = x, y, z), 
normalized to the volume integrated total magnetic energy strength, 
((-B'^)). These quantities are separately averaged. We omit data from 
the smallest domain (FT0.5) for the purposes of clarity The black 
curve is from FT2, red is FT4, blue is FT8, and orange is Y16. As 
with Figure 4, the vertical structure of these quantities are consistent 
between all calculations: the energy density in toroidal fields domi- 
nates throughout the domain. In the region \z\ < 2H, the ratio of 
energy densities is 4 : 40 : 1 (x : j/ : z), while in the region 1^1 > 3H, 
the energy density in vertical fields exceeds that in radial fields, be- 
coming comparable to (but still smaller than) the energy density in 
toroidal fields. 



(e.g. Blandford & Pajme 1982). Repeating run Y16 with an ex- 
tended vertical domain (i.e. a domain size 16H x 32H x 16H) 
could therefore yield insight into the launching of a magnet- 
ically driven wind from a turbulent accretion disc. The com- 
putational cost of such a simulation means that it is beyond 
the scope of the present work. 

4 MESOSCALE SPATIAL CORRELATIONS 

4.1 Magnetic Field 

To examine the structure of the turbulence in the (x, y) plane, 
we employ the autocorrelation function (ACF) (see Guan 
et al. 2009), defined as. 



ACF(/(Aa;)) 



j f{t,x)f{t,x + /^x)d?x 
!f(t,xYd^x 



(8) 



where / is the fluid quantity of interest, and the brackets de- 
note a time-average. Note that we have defined the ACF to be 
normalized by its maximum value (at Ax = Ay — Az = 0), 
and unlike Guan et al. (2009), we do not subtract off any 
mean quantities before calculating the ACF. Furthermore, for 
vector quantities, such as the magnetic field, we take the ACF 
of each component and then sum them as was done in Guan 
et al. (2009); ACF(B) = ACF(_B^) -I- ACF(Bj;) + ACF(BJ. 
The time average is done from orbit 50 to 125 in all cases. 



© 2004 HAS, MNRAS 000, 1-18 



Accretion Disc Turbulence 7 



ACF(B) 



0.0 



0.17 



0.33 



0.50 



0.67 



0.83 



1.0 




5 



-5 


1 1 1 


\ 


1 1 1 








Figure 6. The ACF of the magnetic field, as defined by equation (8) and calculated for \z\ < 2H, for FT2 (left), FT4 (middle), and FT8 (right) 
in the Az = plane. The tilted centroid structure of the magnetic field is the same size and shape in all three domains. 



ACF(B) 




Figure 7. The ACF of the magnetic field, as defined by equation (8) 
and calculated for \z\ < 2H, for the smallest domain (FT0.5, left) 
and the next largest domain (FT2, right) in the Az = plane. While 
the tilted centroid structure of the magnetic field appears to be mostly 
contained within the larger box, it is most definitely not contained in 
the smaller domain. 



Because of the change in vertical structure at |2| ~ 2H, 
we have created two sets of ACFs, one for which we restricted 
the calculation to 1^1 < 2H, the other of which is for [z] > 
2H. This will allow us to probe what has been referred to as 
the turbulent "disc region" {\z\ < 2H) separately from the 
"coronal region" Qz\ > 2H). First focusing on the disc region. 
Fig. 6 shows the ACF of the magnetic field at Az — for the 
runs FT2, FT4, and FT8. To account for the large change in 
scale, we plot the same quantity for the two smallest domains, 
FTO.Sand FT2, in Fig. 7. 

As in the vertical structure plots, there is a striking dif- 
ference between the smallest domain and the larger ones. 
Although for the larger domains, the ACF has a small but 
nonzero value over extended regions (a point that we shall 
return to shortly), it is obvious that the tilted and elongated 
centroid of the ACF is well contained within the domain. This 
is not the case for FT0.5; indeed, it appears that the correla- 
tion length of B extends beyond the size of the region, and 
that the peak of the centroid is not too much larger than the 
value of the ACF away from the centroid. Since Omag is pro- 
portional to the "tilt angle" away from the y axis (Guan et al. 
2009), the low average value of Omag along with the extended 
ACF imply that transport of angular momentum by magnetic 
fields is reliant on the action of strong coherent toroidal mag- 
netic fields, rather than small-scale turbulent structures. 

Beyond the smallest domain, there generally appears to 
be two components to the ACF, one of which is strongly lo- 
calized and tilted and which does not change significantly 
with domain size. The other component is more uniformly 
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Figure 8. The ACF defined by equation (8) for tiie magnetic field 
taken along the major (solid lines) and minor (dashed lines) axes 
of the centroid. The top plot is ACF(_B) calculated for 1^1 < 2H, 
whereas the bottom plot is the same quantity but for |z| > 2H. We 
omit data from the smallest domain (FT0.5) for the purposes of clar- 
ity. The black curve is from FT2, red is FT4, blue is FT8, and orange 
is Y16. The magnetic field in the coronal regions has a larger cor- 
relation length than that within 2H of the mid-plane. In all cases, 
the magnetic field have a strongly localized component as well as a 
non-zero extended component that fills the simulation volume. 



distributed in {Ax, Ay) space and represents the background 
magnetic field. This component fills the entire domain, and in 
that sense cannot converge with increasing domain size. That 
is, as one makes the shearing box larger and larger, there is 
magnetic field that is always volume-filling, and thus there is 
always a component to the field structure that is correlated 
on the largest scales of the box. 

To demonstrate this more quantitatively. Fig. 8 shows ID 
slices of the ACF along the major and minor axes of the tilted 
structure. These plots include the largest box (Y16), but ex- 
clude the smallest domain size. The lower plot corresponds 
to ACF(B) calculated in the coronal region The ACF calcu- 
lation in this region was done for z > 2H and z < ~2H 
separately, and the resulting ACFs were averaged together. 
These plots again demonstrate that there are two components 
to the ACF, one that is concentrated and tilted, the other of 
which is uniform and extended. The tilted correlation struc- 
tures are significantly larger in the coronal region compared 
to the disc region; the magnetic field has a larger correla- 



Figure 9. Time evolution of the volume averaged toroidal field (for 
all x,y and for \z\ < 2H) in code units. FT2 is represented by the 
black curve, FT4 is the red curve, FT8 is the blue curve, and Y16 
is the orange curve. All simulations show the oscillation of the mean 
toroidal field, but the oscillation amplitude is lower in FT8. 



tion length in the coronal region. Furthermore, the value of 
the ACF for the extended component is larger than the cor- 
responding component in the disc region. These observations 
suggest that the mean background field becomes a greater 
fraction of the magnetic energy in the coronal region (this is 
consistent with the two-dimensional correlation functions of 
Guan & Gammie 2011) and that the MRI modes that do op- 
erate in these regions produce fluctuations of relatively large 
length scale. This latter point is likely a result of the charac- 
teristic MRI wavelength becoming larger as the density drops 
(i.e., the Alfvcn speed increases) in the corona. 

Beyond comparing the coronal vs. disc regions, it is worth 
noting the convergence properties of the ACF as the domain 
size is increased. In both the coronal and disc regions, the 
inner centroid appears to be convergent with domain size. 
While the different ACFs at large IS./H do not lie on top of 
each other, they are reasonably consistent with each other, 
with the exception of FT8 (the blue curve). In both the corona 
and disc regions, the extended ACF of FT8 appears signifi- 
cantly lower than the other domain sizes. We have examined 
the evolution of a volume-averaged {By) for these calcula- 
tions (shown in Fig. 9); the FT8 run shows smaller amplitude 
(though still regular) variations in {By), resulting in a smaller 
background field component relative to small scale turbulent 
structures. The smaller background field relative to the tur- 
bulent fluctuations would certainly account for the lower ex- 
tended ACF for FT8. We are not entirely sure why the back- 
ground field is weaker in FT8 in the first place, though some 
statistical variation is to be expected since the evolution of 
{By) is modulated on rather long timescales (Simon et al. 
2011). 

To summarize, the horizontal structure of the magnetic 
field consists of an inner centroid that is tilted with respect 
to the azimuthal direction, as well as an extended compo- 
nent, likely due to the background magnetic field. The gen- 
eral shape of the ACFs converge, though as domain size is 
increased, there will always be structure at the largest scales 
of the box. These results have consequences for the physics of 
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Figure 10. Space-time diagram of fractional gas density fluctuations, <5pfrac (^s defined in the text) in the (t, x) plane. Note the different scales 
on both axes in each plot. For FT4 (upper right), and FT8 (lower left), there exists a kxLx/{2n) = 1 mode zonal flow. The zonal flow in Y16 
(lower right) has a structure more complex than a simple kj:Lx/('2n) = 1 mode. 



angular momentum transport in magnetized accretion discs, 
which we discuss in detail in §4.3. 



4.2 Zonal Flows 

We now turn to the structure and evolution of the gas density 
in the simulated discs. We are particularly interested in the ex- 
istence and properties of zonal flows, long lived axisymmetric 
perturbations to the gas density that result from a geostrophic 
balance between pressure and Coriolis forces. Johansen et al. 
(2009) identified such flows in simulations of MRI disc turbu- 
lence, and suggested that they occur as a result of an inverse 
cascade in magnetic energy, which induces regions of super- 
and sub-Keplerian velocities. In broad terms, zonal flows are 
of physical interest because their existence within protoplane- 
tary discs would provide a way to trap solid particles that are 
partially coupled to the gas disc through aerodynamic forces. 
Here, following prior work by Johansen et al. (2009) and 
Yang et al. (2011), we study whether the size and scale of 
the zonal flows converges with domain size. 

Fig. 10 shows space-time diagrams in the (i, x) plane of 
the fractional perturbation 5pfrac to the density. In calculating 
Sptiac, we average the density over y and z, and then subtract 
off and normalize by the volume averaged p at each time. 
Note the change in both horizontal and vertical scale in each 
panel of the figure. For FT2 (upper left), there is essentially 
no indication of a zonal flow, but for the larger boxes, there 



are large scale long-lived features in p. For FT4 and FT8, 
these features always fill the largest radial scale in the box 
(i.e., a kxLx/{2TT) = 1 mode). The same diagram for the 
largest box (lower right), however, shows a structure more 
complex than a simple k^L^ /{2tv) = 1 mode. 

To quantify the density structure further, we examine the 
ACF of p in the (Ax, Ay) plane (the same averaging is done 
here as for ACF(i3)). For simplicity, we only calculate the 
ACF over l^] < 2H. The result for FT2, FT4, and FT8 is 
shown in Fig. 11, and the ACFs for FT8 and Y16 are shown 
in Fig. 12. Note the range in ACF values from the color bar; 
in general, the density is very nearly uniform in the horizon- 
tal plane. However, in all domain sizes, there is a tilted and 
concentrated structure resembling the magnetic field ACF. 
This particular structure was also seen in the calculations 
of Guan et al. (2009). As domain size is increased, an ax- 
isymmetric structure appears. This structure, which appears 
to be more or less converged between the two largest runs, 
is due to the presence of zonal flows. We have also exam- 
ined a one-dimensional slice along Ay = 0, shown in Fig. 13. 
The evidence for convergence is less persuasive in this one- 
dimensional slice than in Fig. 12, and the most that we can 
conclude is that there are tentative signs that the outer radial 
scale of the zonal flow is Ax ~ 6H (implying a zonal flow 
wavelength of ~ 12H) and is contained within the largest 
spatial domain modeled. Running an even larger domain sim- 
ulation could potentially nail down this radial scale. However, 
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Figure 11. The ACF of the gas density as defined by equation (8) for FT2 (left), FT4 (middle), and FT8 (right) in the Az = plane. In all three 
ACFs, there is a tilted centroid structure. However, as domain size is increased, a larger scale axisymmetric structure appears. 



such large domains become prohibitively computationally ex- 
pensive given our current resources. 

Finally, to explore whether the existence of zonal flows 
is sensitive to the initial conditions, we ran two additional 
simulations. They are the equivalent of FT4 but initialized 
with two flux tubes in one case, and with a uniform toroidal 
field (at constant /?) in the second case. In both cases, a 
kxLj:/{2Tv) = 1 mode zonal flow appears around 10-20 or- 
bits into the simulations. 



4.3 Locality of Angular Momentum Transport 

In the preceding sections, we have demonstrated that there 
are mesoscale correlations within the turbulent disc (\z\ < 
2H) in both the magnetic field and the density. The existence 
of such structures can have important consequences for the 
physics of angular momentum transport in magnetized ac- 
cretion discs. For example, if the mesoscale correlations in 
the magnetic field within the turbulent disc discussed in §4.1 
are actively generated from the magnetorotational instabil- 
ity, then these fields can also transport angular momentum 
through the Maxwell stress term in the angular momentum 
conservation equation (Balbus & Papaloizou 1999). Similarly, 
the zonal flows discussed in §4.2 are associated with sub- and 
super-Keplerian velocities which form due to the action of 
magnetic forces (Johansen et al. 2009); such flows can then 
be associated with transport of angular momentum as they 
drift through the disc. In both of these cases, we would ex- 
pect to measure angular momentum transport (through the 
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Figure 12. The ACF of the gas density as defined by equation (8) 
for FT8 (left) and Y16 (right) in the Az = plane. The density 
structure is consistent between these two domain sizes, including the 
large scale axisymmetric component. 



Maxwell and Reynolds stress, respectively) on scales 2> H. 
Such a result would undercut the assumption that angular 
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Figure 13. The ACF of the gas density as defined by equation (8) 
versus Ax along Ay = and Az = 0. The colors denote the vari- 
ous runs; black is FT2, red is FT4, blue is FT8, and orange is Y16. 
All ACFs show a central component that flattens out at ~ 0.3H. For 
boxes larger than FT2, there is a second dip corresponding to the 
presence of a zonal flow. Finally, there is hint of this second dip flat- 
tening out for Y16, thus providing tentative evidence for an outer 
scale to the zonal flow that is contained within the simulation do- 



momentum transport in magnetized discs can be treated as a 
purely local process (Shakura & Sunyaev 1973). 

We can test these ideas by examining the ACF for the 
Maxwell and Reynolds stress. If angular momentum trans- 
port within the disc is determined locally, we would expect 
the ACF of the stress to decay rapidly on scales > H; cor- 
relations on scales 3> H would indicate a non-local compo- 
nent to angular momentum transport (Gammie 1998). Fig- 
ure 14 plots one-dimensional slices along the major and mi- 
nor axes of the tilted centroid component of the ACF for 
the Maxwell stress, ACF{—BxBy) (top panel), and Reynolds 
stress, ACF{pvx5vy) (bottom panel), calculated over the re- 
gion \z\ < 2H, i.e. within the "turbulent disc". The data of 
this figure demonstrate that while the ACFs are strongly con- 
centrated within < H, there is an extended component, cor- 
related at ~ 20% for the Maxwell stress and ~ 8% for the 
Rejmolds stress, that fills the simulation domain. Essentially, 
these results imply that there is a net background Maxwell 
and Reynolds stress resulting from mesoscale structures in the 
magnetic field and density, on top of which turbulent stress 
fluctuations are imposed. 

To verify this result, we (for all runs except FT0.5) exam- 
ined the time-averaged (from orbit 50 onward) vertical pro- 
file of the quantity 



1.00 



(9) 



{-BxBy) 

which represents the fraction of the average Maxwell stress 
that resides in the large scale background magnetic field. We 
find that this ratio is quite small near the disk mid-plane but 
increases dramatically away from the mid-plane, leveling out 
to ~ 0.3—0.4 for \z\ > 2H. The conclusion from these calcula- 
tions is clear: while angular momentum transport within mag- 
netized accretion discs is highly localized, there is, in addi- 
tion, a significant non-local component, the existence of which 
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Figure 14. As in the top panel of Figure 8 for the ACF of the Maxwell 
(top panel) and Reynolds (bottom panel) stress. While the ACFs are 
strongly concentrated within < H of the centroid, there is an ex- 
tended component, correlated at ~ 20% for the Maxwell stress and 
~ 8% for the Reynolds stress, that fills the simulation domain in all 
cases. 



implies that models founded upon a purely local description 
of angular momentum transport are incomplete. 

The origin of these mesoscale correlations in the Maxwell 
stress can be understood, at least in part, by considering the 
saturation characteristics given in Table 2. In all of the sim- 
ulations, the turbulent disc is dominated by strong toroidal 
magnetic fields (making up ~ 80% of the magnetic energy in 
simulation Y16). The toroidal field MRI is well-resolved ac- 
cording to the criteria of Sano et al. (2004), which specifies 
the number of zones per characteristic wavelength of the in- 
stability. In principle, all modes with wavelengths longer than 
this characteristic wavelength are unstable. This is in contrast 
to the case of the vertical field MRI, where wavelengths longer 
than that of the pressure scale height are expected to be sta- 
bilized. Thus, it is not unexpected that long wavelength struc- 
tures form in the turbulence from the toroidal field MRI. We 
discuss this further in Section 7. 



5 TEMPORAL VARIABILITY 

We now return to the behaviour of the temporal variability in 
domains of different size. From Fig. 2, it is clear that there is 
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Figure 15. Space-time diagram in (t, 2) for tiie liorizontally averaged By in FT2 (upper left), FT4 (upper right), FT8 (lower left), and Y16 
(lower right) . The so-called "butterfly diagram" present in vertically stratified MRI-driven discs is apparent in all four domain sizes, and becomes 
more regular as the domain size is increased. 



a marked decrease in temporal variability as the domain size 
is increased. This is generically true for a variety of quantities. 
Figure 15, for example, shows the (t, z) space-time diagram 
of By for the four largest domain sizes. The dominant feature 
in this diagram is the ~10 orbital period flipping of By that 
has been observed in many previous calculations (e.g., Bran- 
denburg et al. 1995; Davis et al. 2010; Simon et al. 2011; 
Guan & Gammie 2011). However, superimposed upon this is 
a stochastic component to By, which becomes smaller as the 
domain size is increased. 

To better quantify this behavior, we define a diagnostic, 
which represents the fractional power in temporal fluctua- 
tions. 



Var[{/)] 
Avg[(/)]2 



(10) 



Here, "Var" refers to the temporal variance of the quantity 
of interest, and "Avg" is the temporal average of the quan- 
tity. Figure 16 shows this quantity versus box size for both 
the total stress (squares) and the magnetic energy (aster- 
isks), calculated from orbit 50 onward within the "turbulent" 
disc i\z\ < 2H). The data of this figure demonstrate that 
there is a clear decrease in the variability with increasing do- 
main size. Specifically, as we increase the domain size from 
2H xAH X 8H (FT2) to 8H x 16H x 8H (FT8) (a factor 16 
increase in volume), we find that the fractional variability de- 
creases by approximately a factor 10. Increasing the volume 
by a further factor of four (i.e. from simulation FT8 to Y16) 



1.000 


m 










■ 


0.100 


□ 


n 


g 






■ 


0.010 


^ 






X 

□ 




^ 


0.001 















10 



L,/H 



Figure 16. Variability of the volume-averaged stress (squares) and 
magnetic energy (asterisks), as defined by Pvar (equation 10) and 
averaged over all x and y and for 1^1 < 2H, versus box size. As the 
domain size is increased, the average temporal variability decreases. 



results in only a (comparatively) small decrease in fractional 
variability, suggesting that we are approaching convergence 
at this largest domain size. 

One might expect that the volume averaged variability 
will decrease as domain size is increased; as we average over 
increasing numbers of (largely) uncorrelated subvolumes, the 
fractional variability of the system will decrease. If this sta- 
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Figure 17. Variability of tlie stress (red squares) and magnetic energy (black asterisks) as defined by Pvar (equation 12) versus box size. This is 
the same quantity as is plotted in Fig. 16. However, the total domain is broken into horizontal subdomains of size 2H x 4H, after which Pvar is 
calculated for that subdomain. All values of Pvar are then averaged together and plotted against the box size. The error bars denote one standard 
deviation about the average over subdomains. As labelled, each panel corresponds to a different vertical region in the shearing box. The primary 
result of this figure is that near the mid-plane, temporal variability decreases uniformly throughout the grid as domain size is increased, whereas 
for 1^1 > 2H, this variability is constant with domain size. 



tistical effect is solely responsible for the observed decrease 
in variability, we would expect the variability associated with 
a small domain within the simulation volume to remain un- 
changed as the volume is increased. To test this notion, we 
divide each shearing box into horizontal subdomains of size 
LxXLy — 2Hx4:H (i.e. the domain size of simulation FT2 on 
the X — y plane), calculate Pvar in each subdomain, and then 
average the Pvar. This quantity is plotted in Fig. 17 for a va- 
riety of vertical regions in each shearing box, excluding the 
smallest, FT0.5. The error bars denote one standard devia- 
tion about the average Pvar value. Thus, the error bars do 
not represent temporal variability but instead quantifies the 
spread in Pvar values between the subdomains. 

In the corona \z\ > 2H, the variability of the subdomains 
remains constant as domain size increases, but in the turbu- 
lent disc (|z| < 2H), the variability of the subdomains de- 
creases with increasing domain size. As the shearing boxes are 
increased in size, the physics of the turbulence changes such 
that variability within the turbulent disc (1^1 < 2H) decreases 
everywhere uniformly, while the variability in the coronal re- 
gion remains constant. 



COMPARING THE MESOSCALE AND GLOBAL 
REGIMES 



Measured in units of H, our largest shearing boxes span ra- 
dial extents that are not so different from those modeled in 
some global calculations. It is therefore possible to directly 
compare some of our results to those obtained from stud- 
ies of disc turbulence in the global regime. We restrict our 
comparison to recent global simulations which attain com- 
parable resolutions to the 32 zones per H of the local runs 
(such as Beckwith et al. 2011; Sorathia et al. 2011). The 
most secure comparisons are probably with the simulations 
of Sorathia et al. (2011); these have orbital advection imple- 
mented, which leads to a significant reduction in numerical 
diffusion. We note that quantitative comparisons are not al- 
ways possible: there are significant and unavoidable differ- 
ences in the setup and analysis between the local and global 
runs. We also note that the scale height defined in Beckwith 
et al. (2011) and Flock et al. (2011a) is a factor of \/2 smaller 
than the scale height defined here. While this is an order unity 
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Figure 18. Time evolution of thie volume averaged toroidal field in 
code units (for all x,y and for |2| < 2H). The solid black curve is 
FT8, the dashed curve is Y16, and the red curve is global simula- 
tion data from Beckwith et al. (2011) (their Fig. 3), renormalized to 
match the vertical scale of the shearing box data and plotted against 
the time unit of the local simulations. The early evolution of the net 
toroidal flux is comparable between the global and local simulations; 
the initial net flux is rapidly expelled after which the flux oscillates 
about zero. The original global simulation data from Beckwith et al. 
(2011) is for 20 orbits at 15 Schwarzschild radii. The agreement be- 
tween the curves in the figure suggest that the dynamo present in 
the simulation of Beckwith et al. (2011) has a characteristic radius 
of 10 Schwarzschild radii, at which the duration of each cycle of the 
dynamo is 10 orbits and agrees with local simulations. 



effect, it should be borne in mind when making direct com- 
parisons between that work and the work presented here. 

To start, we consider the standard a and Qmag parame- 
ters. In all of our calculations, we find that a ~ 0.02 — 0.03, 
not significantly different than the same parameter calculated 
(albeit calculated differently) in global simulations. Figure 5 
of Beckwith et al. (2011) displays a versus time. While there 
is a long term decrease in the a value, their values range 
from 0.02 to 0.06. The results of Flock et al. (2011a) give 
a ~ 0.005 — 0.02, somewhat lower than our values and the 
Beckwith et al. (2011) values. Sorathia et al. (2011) calcu- 
lated a for various initial field geometries; they find that a 
ranges from ~ 0.01 to ~ 0.08 depending on the field geome- 
try and resolution. 

Again considering Fig. 5 of Beckwith et al. (2011), the 
dashed lines in the bottom panel show Omag as a function 
of time, which is relatively constant after the initial growth 



of the MRI. The values of . 



range between 0.3 and 0.4, 



generally consistent with all of our runs (except for FT0.5 
as elaborated upon above), for which Qmag ~ 0.4. The high- 
est resolved global simulations in Hawley et al. (2011) have 
Omag ~ 0.2—0.4, again generally consistent with our shearing 
box results (though, somewhat lower). This quantity can be 
related to the tilt angle of the magnetic correlation function 
(Guan et al. 2009; Sorathia et al. 2011), which is what So- 
rathia et al. (2011) focused on calculating. Our typical value 
of Qmag ~ 0.4 corresponds to 9tat ~ 12° in general agreement 
with the Sorathia et al. (2011) calculations. The global calcu- 
lations of Flock et al. (2011a) return Omt ~ 8 - 9°, slightly 
lower than our value as well as that of Sorathia et al. (2011). 



Figure 19. Volume averaged stress normalized by the volume aver- 
aged gas pressure versus time for FT8 (black solid line), and Y16 
(dashed line). This average was done over all x and y and for 
\z\ < H. The red curve is the equivalent normalized stress from 
Beckwith et al. (2011). The horizontal scale is chosen to match that 
of Fig. 18, and as such, the global data is more representative of a 
radius of 10 Schwarzschild radii, as explained in the text and Fig. 18. 
Similar variability levels are observed in both local simulations and 
global simulations. 



Other behaviors are best compared in a more qualitative 
manner. The "butterfly diagram" appears to be a robust fea- 
ture of vertically stratified MRI-driven turbulence. It has been 
seen in several global calculations (e.g., O'Neill et al. 2011; 
Beckwith et al. 2011; Flock et al. 2011b,a), and the period of 
the azimuthal field sign flipping is approximately 10 local or- 
bits with power shared radially (O'Neill et al. 2011). Looking 
again at the time histories of the pressure-normalized total 
stress in global simulations, the temporal variability in the 
global simulations appear to be quite small (Beckwith et al. 
2011; Flock etal. 2011a), similar to the largest shearing boxes 
presented here. A more explicit comparison can be made by 
examining the history of the largest shearing boxes presented 
here over the duration of a well-resolved global simulation. 
To do so, we note that Fig. 3 of Beclavith et al. (2011) plots 
the time-history of the toroidal magnetic flux within the sim- 
ulation domain, which exhibits transient behaviour, followed 
by two complete dynamo cycles. In Fig. 18, we combine this 
global simulation data with the history of simulation FT8 and 
Y16. Note that the local data is plotted over 40 orbital periods, 
whereas the data presented by Beckwith et al. (2011) is for 
20 orbits at 15 Schwarzschild radii. This would suggest that 
the dynamo present in the simulation presented in Beckwith 
et al. (2011) has a characteristic radius of 10 Schwarzschild 
radii, at which the duration of each cycle of the dynamo is 
10 orbits. The figure shows a comparable evolution between 
the local and global simulations; transient behavior (corre- 
sponding to expulsion of toroidal flux from the simulation do- 
main) lasts for ~ 15 orbits, after which the MRI dynamo pro- 
duces toroidal magnetic fields of alternating sign on a 10 orbit 
timescale. Figure 19 shows the stress evolution comparison 
over this same time period for FT8 and Y16 averaged over 
all X and y and ior \z\ < H. A more direct comparison will be 
necessary in the future (e.g., calculating the Pvar parameter 
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in a shearing box sized region of a global simulation), but this 
basic comparison suggests that our largest boxes are reaching 
the temporal variability level seen in global simulations. 

The vertical structure of the turbulence, namely the mag- 
netic field, is another point of comparison. Our calculations 
show a roughly flat distribution of the gas pressure normal- 
ized stress at ~ 0.01 for \z\ < 2H, after which the stress drops 
to ~ 10^'' (see Fig. 4). The same figure shows that /9 peaks 
near 40 and then drops below unity outside of |z| = 2H. 
In Beckwith et al. (2011), the normalized stress peaks near 
~ 0.01 and /3 ~ 30 at its peak, generally consistent with our 
results. However, the behavior away from the mid-plane dif- 
fers dramatically from what is seen in our local calculations. 
In particular, the stress does not drop nearly as rapidly out- 
side of l^l = 2H, and /? never drops below unity. This dif- 
ference in behaviors away from the mid-plane could be due 
to the use of periodic boundary conditions in the global sim- 
ulations of Beckwith et al. (2011). The general structure of 
/3 in the global BO model of Flock et al. (2011b), which has 
vertical outflow boundary conditions, (their Fig. 11) is similar 
to what we observe in our local calculations except that ev- 
erything appears to be scaled up by a factor of ~ 10 in Flock 
et al. (2011b). That is, their peak j3 value is ~ 500, and at 
l^l « 2H, P ~ 10. The reason for this large difference is not 
at all clear. 

All of these comparisons, with the exception of the ver- 
tical magnetic structure, yield good agreement between local 
and global calculations, considering the spread in resolution 
and measurement techniques. However, the real question we 
have set out to answer in this paper is whether or not the 
mesoscale appears the same between these two calculations; 
this has yet to be addressed. The most relevant measure for 
answering this question is the correlation function of various 
turbulent structures in the horizontal plane. 

The analysis of Beckwith et al. (2011) and Flock et al. 
(2011a), shows a tilted, highly localized component to turbu- 
lent fluctuations, in agreement with the results found in this 
work. In their analysis, Beckwith et al. (2011) subtract off 
large scale variations (due to radial and vertical gradients) to 
only focus on the central, tilted component. They then cal- 
culate a correlation length along the major and minor axes 
of the tilted structure by determining the length at which the 
correlation function drops by a factor of e~^ (see their Ta- 
ble 1). Flock et al. (2011a) calculates the same lengths but 
by using the half width at half maximum of the correlation 
function. The correlation length values between these two 
papers generally agree, though it should be noted that Beck- 
with et al. (2011) consider the ACF of the magnetic energy, 
whereas Flock et al. (2011a) calculate the ACF for B as was 
done here and in Guan et al. (2009). 

In our calculation of the ACF, we do not subtract off any 
background component. However, we can consider the corre- 
lation length along a given axis to be the point at which the 
ACF begins to flatten out. Beyond this length, the ACF is flat 
because the structure is relatively uniform; thus, the "knee" in 
the curve represents the outermost scale of the localized com- 
ponent to the ACF. From Fig. 8, it appears that the correlation 

-0.3 -0.4//. 



lengths for B axe. \b,, 



4 — 5H and As.i 



Viewing the equivalent plot for B and p gives Afl2 j^j^j ~ 
2 - 3H, As2,„,i„ ~ 0.2H, and Ap,min ~ 0.3 - OAH. We did 
not calculate a major axis correlation length for p because the 



ACF never flattened out along the major axis. This is likely 
related to the presence of the axisymmetric zonal flows enter- 
ing into the ACF calculation, though it is still not clear why 
the ACF does not flatten out at some point along the major 
axis. In any case, it would seem that there is rough agree- 
ment between the correlation length of the tilted, localized 
ACF component in global and local calculations. 



The general agreement between these correlation 
lengths is very encouraging as it shows that MRI-driven tur- 
bulence has the same small scale structure in both local and 
global simulations. However, a primary result of the work 
here is that in addition to the tilted component of the mag- 
netic field structure, there is also a non-negligible, volume 
filling background component. Again, since the large scale 
component was subtracted off in the analysis of Beckwith 
et al. (2011), a direct comparison with our results cannot 
be made very easily. However, some information can be ex- 
tracted by examining their Fig. 4, which shows the structure 
of the toroidal magnetic field on the azimuthal plane. The 
figure shows that this field has structure in both the vertical 
and radial dimensions; there appear to be magnetic field bun- 
dles that evolve throughout the disc. Through a very rough 
examination of this figure, it appears that the radial size of 
these bundles is ~ 5rs, where Vs is the Schwarzschild ra- 
dius. In those calculations, H/R = 0.07 so that at _R = lOr^, 
H = 0.7rs. Thus, the tjqsical size of the magnetic field bundles 
in that figure is estimated to be ~ 7H where (as noted above) 
this H is that defined in Beckwith et al. (2011), a factor of \/2 
smaller than the H that we use. This scale translates to ~ 5H 
using H defined as we have in this paper, which should defi- 
nitely be captured in our largest shearing box, Y16. However, 
the correlation function of B clearly indicates that there is a 
volume-filling background component to the magnetic field; 
there is no structure associated with ~ 5H. Granted, this com- 
parison is murky at best, and perhaps even larger shearing 
boxes would show the presence of such radial structure. How- 
ever, these results suggest that there may be a fundamental 
difference in radial magnetic structure at the mesoscale be- 
tween local and global simulations, perhaps a result of the 
inclusion of curvature terms in global simulations. 



The presence of density zonal flows in global calculations 
done to-date is questionable at best. To our knowledge, the 
best evidence for such flows exist in the global calculations of 
Lyra et al. (2008); see their Fig. 10, which shows reasonably 
long-lived radial variations in the gas density. However, even 
here, the flows are only weakly detected and not discussed 
by the authors in any detail. It would thus seem that these 
flows are very difficult to detect in global calculations. This 
difficulty may result from larger scale density gradients domi- 
nating the weaker zonal flow variations. While our work here 
has shown tentative evidence for a converged outer scale to 
these flows in shearing boxes, it is still not clear whether or 
not such flows are an artifact of the shearing box approxima- 
tion that tend to exist on the largest radial scales. Thus, we 
believe that the detection of these structures in global simula- 
tions should be a very important priority for future work, one 
that we are currently pursuing. 
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Figure 20. Fraction of the total magnetic energy as a function of 
isotropic wavenumber k. In calculating these curves, we averaged 
the three-dimensional Fourier transform of magnetic energy over 
spherical shells of constant k and then averaged the resulting one- 
dimensional power spectra in time, from orbit 50 to 125. We omit 
data from the smallest domain (FT0.5) for the purposes of clarity. 
The black curve is from FT2, red is FT4, blue is FT8, and orange is 
Y16. As domain size is increased, the fraction of the total magnetic 
energy at the smallest scales decreases while the largest fraction of 
magnetic energy resides at the largest scales. There appears to be 
convergence between the two largest domain sizes. 



7 DISCUSSION, SUMMARY, AND CONCLUSION 

We have examined the convergence properties of MRI-driven 
MHD turbulence with increasing domain size, i.e. as we tran- 
sition from the "local" (if ~ i < 71) to "mesoscale" (H < 
L <^ R) regime. We have found that volume-averaged quan- 
tities (such as the Shakura & Sunyaev 1973, a-parameter and 
the ratio of the Maxwell stress to magnetic energy, Qmag) 
converge rapidly as we transition between these regimes. We 
have also studied the properties of the magnetic field autocor- 
relation function (ACF) as we move between these regimes, 
which reveals a two component structure to the magnetic 
field. One is highly localized and tilted with respect to the 
azimuthal direction, consistent with previous findings (e.g., 
Guan et al. 2009). The other component is extended to the 
largest scales in the domain, and is likely associated with the 
(predominantly toroidal) background magnetic field. Further- 
more, we have identified the zonal flows of Johansen et al. 
(2009) in our domain sizes of scale AH x d,H x 8H and 
larger These flows always fill the largest radial scale in the 
box except for our largest simulation. We find tentative evi- 
dence that the radial wavelength of these flows converge at 
~ 12H. 

These results have important implications for the physics 
of angular momentum transport in magnetized accretion 
discs. While the ACFs of the Maxwell and Reynolds stress are 
strongly concentrated within < H, both of these quantities 
are also correlated on scales ^ H at the ~ 20% and ~ 8% 
level respectively. If angular momentum transport in magne- 
tized accretion discs is a purely local process, then the ACF of 
these quantities should drop to zero on scales ~ H (Gammie 
1998). That the ACFs of the accretion stress remain finite at 
scale ^ H undercuts the assumption that angular momen- 



tum transport in magnetized discs can be treated as a purely 
local process (Shakura & Sunyaev 1973; Pringle 1981). 

We have also found that the transition to the mesoscale 
limit is associated with a decrease in the temporal variabil- 
ity of the turbulent disc, an effect that is not solely due to 
volume-averaging over uncorrelated subdomains. It is per- 
haps not surprising that the small scale physics of MRI-driven 
turbulence changes with increasing box size; as larger mag- 
netic structures are permitted in the box, large scales can po- 
tentially interact with smaller scales and change the variabil- 
ity properties at these scales. While a spectral energy transfer 
analysis (Fromang & Papaloizou 2007; Fromang et al. 2007; 
Simon & Hawley 2009) would be useful in understanding 
how exactly different scales interact, we can attain a basic 
understanding by considering the time-averaged power spec- 
tra of the magnetic energy, as shown in Fig. 20. In calculating 
this power spectra, we took a full three-dimensional Fourier 
transform, averaged the result over spherical shells of con- 
stant k, and then time averaged the resulting spectra from 
orbit 50-125. The figure shows the fraction of the total mag- 
netic energy that exists at a scale k (i.e., integrating each 
curve over k would return unity) . As box size is increased, the 
power spectra shift downwards because there are more spa- 
tial scales over which the magnetic energy can be distributed. 
However, the main point to take away from this figure is that 
as the box size is increased, the fraction of the total mag- 
netic energy at the smallest scales generally decreases while 
the largest fraction of magnetic energy resides at the largest 
scales, with signs of convergence between the two largest do- 
main sizes, (this result is consistent with Fig. 2 of Johansen 
et al. 2009). Since the turbulent energy levels are roughly the 
same between all of these simulations (e.g.. Fig. 3), increas- 
ing the domain size equates to more power being taken from 
the small scales and put into the largest scales. How and why 
exactly this behavior occurs is still an open question. How- 
ever, this result is consistent with our notion that the largest 
scales in the shearing box play a very important role in the 
evolution of the MRI. 

We have also examined the physical properties of shear- 
ing boxes computed using restricted spatial dimensions, e.g. 
0.5H X 2H X 8H in L^x LyX L^. Turbulence arising from the 
MRI in these simulations is characterized by a greater level of 
isotropy (omag = 0.16) than is found in larger shearing box 
domains (those with L^ > 2H x Ly > AH). Furthermore, 
we find anomalous structures in both the ACFs and vertical 
profiles of (e.g.) the magnetic energy. These results suggest 
that simulation domains of this size or smaller could produce 
misleading results. One such use of small domains (neces- 
sary for numerical reasons) is in radiation pressure dominated 
MRI simulations (e.g., Hirose et al. 2009; Blaes et al. 2011). 
However, the ACF of the magnetic field in these simulations 
strongly resembles that of our converged FT2 calculation (Hi- 
rose, private communication), despite the fact that the scale 
height of these radiation pressure dominated simulations is 
not very different from an isothermal disk scale height as 
we have in our calculations (Krolik, private communication) . 
This surprising result indicates that the thermodynamic prop- 
erties of the disk may play a role in determining the detailed 
structure of MRI-driven turbulence. 

We have compared results from the largest domain sizes 
considered here (fiH x l&H x 8H and 16H x 32H x 8H) 
with those obtained from well-resolved global simulations 
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that have been recently presented in the Hterature (e.g., Beck- 
with et al. 2011; Sorathia et al. 2011). We have found broad 
agreement between these two classes of simulations, in terms 
of volume-averaged quantities, the vertical structure of and 
correlations lengths within the turbulence. A point of con- 
tention however, is the largest spatial scale on which (e.g.) 
the magnetic field is correlated. The mesoscale simulations 
presented here suggest that correlations exist on radial scales 
up to the maximum considered (16 H). Crude analysis of the 
results of Beckwith et al. (2011) suggests, however, a maxi- 
mum radial correlation length of ~ 5H. This result may indi- 
cate that the maximum radial correlation length of the mag- 
netic field is determined by global (e.g., curvature), rather 
than local effects. Furthermore, if mesoscale structures in 
the magnetic field and accretion stress originate in the ac- 
tion of the toroidal field MRI, then global simulations that 
probe up to m = 1 in azimuth, utilizing resolutions equiva- 
lent to the simulations presented here, will be necessary to 
determine the outer scale of these quantities. Since there are 
non-negligible levels of angular momentum transport at the 
largest azimuthal scales, a first-principles global treatment of 
the system is necessary to understand the complete physics of 
angular momentum transport. 

Overall, our results suggest several avenues of future 
work. If, as suggested by our results and those of Beckwith 
et al. (2011), a significant amount of angular momentum 
transport takes place on scales 3> H, it will be important 
to determine the length scale on which the work done by 
this large scale angular momentum transport is dissipated. 
Hirose & Turner (2011) suggest that the magnetic fields as- 
sociated with this large scale angular momentum transport 
can rise buoyantly into the disc corona (1^1 > 2H), where 
they can be dissipated. If this is the case, then the local con- 
nection between angular momentum transport and disc heat- 
ing would be broken, further undermining the assumptions 
of local models of angular momentum transport (Balbus & 
Papaloizou 1999). It would then be important to understand 
the mechanisms by which the dissipated heat would be radi- 
ated away by the disc and any resulting changes to the vertical 
structure. 

Equally important will be to identify the physical mech- 
anism that causes the decrease in turbulent variability as we 
transition to the mesoscale regime. As it currently stands, our 
work suggests that the variability of an ionized accretion disc 
is tied to the global structure of the flow. These results, in 
combination with our results regarding the locality of angular 
momentum transport, suggests that attempts to interpret ob- 
servations of accretion variability through local, viscous mod- 
els of angular momentum transport (King et al. 2007; Ingram 
& Done 2011) could yield misleading conclusions regarding 
the physics of magnetized accretion discs. Future global simu- 
lations that incorporate boundary conditions and physics ap- 
propriate to a given astrophysical system are likely to be nec- 
essary to provide accurate interpretations regarding the vari- 
ability of magnetized accretion discs. 

Finally, the detection of density zonal flows in global cal- 
culations will be a very important avenue to pursue. As men- 
tioned above, it will be necessary to determine whether these 
flows are an artifact of the shearing box approximation or are 
in fact present in global simulations as well. If they do exist 
in global calculations, we will need to quantify their radial 
(and azimuthal) scale as well as their amplitude relative to 



the background density. Such quantifications will be impor- 
tant in informing models of planetesimal trapping and planet 
migration. 
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